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ABSTRACT 

We present the results of a study that uses numerical simulations to interpret 
observations of tidally disturbed satellites around the Milky Way. When analysing 
the simulations from the viewpoint of an observer, we find a break in the slope of 
the star count and velocity dispersion profiles in our models at the location where 
unbound stars dominate. We conclude that 'extra-tidal' stars and enhanced velocity 
dispersions observed in the outskirts of Galactic satellites are due to contamination 
by stellar debris from the tidal interaction with the Milky Way. However, a significant 
bound population can exist beyond the break radius and we argue that it should not 
be identified with the tidal radius of the satellite. 

We also develop and test a method for determining the mass loss rate from a Galac- 
tic satellite using its extra-tidal population. We apply this method to observations of 
globular clusters and dwarf spheroidal satellites of the Milky Way, and conclude that 
a significant fraction of both satellite systems are likely be destroyed within the next 
Hubble time. 

Finally, we demonstrate that this mass loss estimate allows us to place some limits 
on the initial mass function (IMF) of stars in a cluster from the radial dependence of 
its present day mass function (PDMF). 

Key words: globular clusters: general - stellar dynamics 



1 INTRODUCTION 



Ian fc Heggie 1997 ; Elson, Hut fc Inagaki 1987). Ideally we 



Accurate estimates of mass loss rates from the satellite 
galaxies and globular clusters in orbit around the Milky Way 
offer the possibility of better understanding the dynamical 
history of the Galaxy. When integrated over the entire pop- 
ulation, such measurements would tell us about the current 
accretion rate onto the Galaxy. This in turn would indi- 
cate whether the Galaxy has grown significantly through the 
gradual process of tidal stripping and disruption of its satel- 
lites. The direct inference of mass loss rates from individual 
objects would place new constraints on detailed analytic or 
numerical dynamical models of these systems. In addition, 
the measurements could be used to identify satellites which 
are likely to have well-populated streams of tidal debris as- 
sociated with them. 

In particular, the formation, evolution and ultimate fate 
of the Galactic globular clusters is a long standing puzzle in 
astrophysics. Dynamically, globular clusters are clean and 
relatively simple systems, yet on close inspection they re- 
veal a wealth of complex behaviour, including relaxation and 
evaporation from internal dynamical interactions and mass 
loss in response to the tidal field of the Milky Way (Mey- 



would like to be able to assess the importance of each of 
these processes for the evolution of individual globular clus- 
ters from observations. 

Of particular interest for mass-segregated, relaxed clus- 
ters is the differential mass loss of stars due to the differ- 
ent relative densities of low and high mass stars at various 
radii, and the consequence of this loss for the evolution of 
the present day mass function (PDMF) away from the initial 
mass function (IMF) . Understanding the form of the IMF of 
the current globulars is relevant for theories of star f orma- 
tion and globular cluster formation ( ^all fc Rees 1977 ), and 
the related issue of how the Galaxy was assembled. The mass 
function of a stellar system dN/dM is defined as the number 
of stars dN in the mass interval (M, M + dM). A common 
functional form chosen to approximate this distribution is 



dN 
'dM 



oc M' 



-(x + l) 



(1) 



where x is called the mass function index. The PDMF x 
of a cluster has been found to be related to the Galacto- 
centric radius Rgc and height above the Galactic plane Z 
of the cluster (Capaccioli, Ortolani fc Piotto 1991). Using 
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simple semi-analytic models Capaccioli, Piotto fc Stiavelli has received r e newed attention (Weinberg 1994 ; Kundic & 



(1993) ihowe^^hat_^h£se_correlato Johnston et al. 1998; Gnedin, Hernquist & 



by evol ring a population of globular clusters which all have 
the same initial IMF in the tidal field of the Milky Way. 
Similarly, in a comparison of several globular clusters Pi- 
otto, Cool & King (1997) found the luminosity function of 
NGC 6397 to be much flatter than M15, M30 and M92, and 
concluded that this could be due to the extreme dynami- 
cal evolu tion implied by its orb it calculated from its proper 
motion ( Dauphole et al. 1996 ). In addition McClure et al. 
(1986) reported a dependence of the PDMF on metallicity, 
and Djorgovski, Piotto & Capaccioli (1993) successfully in- 
corporated this together with the trends in -Rgc and Z into 
a trivariate analysis. In contrast, in a recent study of Pal 
1, Rosenberg et al. (1997) found its mass function to be in- 
consistent with the trends in iicc, Z, and metallicity and 
concluded that this could either be due to evaporation of 
low mass stars, or an intrinsically different IMF. 

One step towards solving the puzzle of just how far 
the PDMF of a cluster differs from its IMF is determining 
how to interpret the signatures of tidal effects in observa- 
tions of globular clusters and other satellites. In the past, 
the limiting radii of Galactic satellites (which are assumed 
to correspond to the tidal radii nide imposed by the Milky 
Way) have been interpreted using King's (1962) tidal radius 
formula 



Ostriker 1998). However, all these studies predict that the 



r-tidc = ^GC 



(■ 



V3Mgc 



(2) 



where Rgc is the distance of the satellite from the Galactic 
centre, rriaat is the mass of the satellite, and Mgc is the mass 
of the Galaxy enclosed within i?GC- For the observed rtido 
this f ormula has been us ed to: (i) find rrisat from the current 
_Rgc (Faber & Lin 1983); (ii) estimate the pericentre of the 



satellite's orbit given the value for rrisat implied by its inter- 
nal dynamics (Oh, Lin & Aarseth 1995; Irwin & Hatzidim- 
itriou 1995 - hereafter IH85); and (iii) com pare t he Galaxy's 
tidal field with the satellite's internal field ( IH95 ) . Such anal- 
yses have been complicated by the detection of 'extra-tidal' 
stars around both dwarf spheroidal galaxies (IH95; Kuhn, 
Smith & Hawley 1996) and globular clusters (Grillmair et 
al. 1995 - hereafter G95; for a summary of current observa- 
tions of both Galactic and extra-galactic globular clusters, 
see Grillmair 1998), and these morphological features have 
been shown to be consistent with those seen in simula t ions of 
tidally disrupting systems ( |0h, Lin fc Aarseth 1995 ; G95). 



Hut & Djorgovski (1992) also estimated destruction 
rates of globular clusters by finding the extent to which their 
distribution of half-mass relaxation times deviates from a 
power law and interpreting this difference in terms of evolu- 
tion of the system. 

An alternative to these very direct interpretations of 
observations is to build models that include some represen- 
tation of all the expected dynamical effects on globular clus- 
ters and inte grate the cluster's e v olution over the lifetime of 
the Galaxy (|Aguilar et al. 1988|; [Capriotti fc Hawley 1996t 



1997 



Gnedin fc Ostriker 1997 ; |Murah fc Weinberg 1997 f Vesperini 



The success of these semi-analytic methods (as with 



those in the previous paragraphs) rests upon the representa- 
tion of the physics involved, and there is some disagreement 
about the level of sophistication required for accurate pre- 
dictions - one example is the theory of tidal shocking which 



Milky Way's globular cluster system is currently undergoing 
significant evolution. 

In this paper, we return to the earlier philosophy of try- 
ing to understand how much we can learn directly from ob- 
servations of tidal signatures. We approach the problem by 
using numerical simulations to assess how successfully such 
features can be interpreted. Our study differs from earlier 
ones in several ways: (i) we include a mass spectrum, and 
begin with mass segregated models; (ii) there is a one-one 
representation of stars in e.g. globular clusters; and (iii) the 
tidal field of a full Milky Way model is included rather than 
idealised as giving rise to a spherical tidal boundary or being 
represented by a one-component potential. In our analysis, 
we explicitly distinguish between the bound and unbound 
stars in a satellite, which allows us to identify the character- 
istics of the extra-tidal population. We use this analysis to 
develop and test methods for quantifying the mass loss rate 
and evolution of the mass function of clusters and satellites 
directly from current observations. 

We present the simulation and analysis methods used 
throughout the paper in §2. The simulations are used to pro- 
vide an overview of the characteristics of evolution in a tidal 
field in §3. We analyse our models from the viewpoint of an 
observer §4. In §5, the interpretations developed in §4 are 
applied to observations and used to measure the rate of de- 
struction of the Milky Way's satellite system. We summarise 
our results in 56. 



2 METHODS 

2.1 General approach 

In our calculations, the satellite is represented as a discrete 
set of particles of different masses, one for each star expected 
in a cluster of the chosen mass and mass function. The par- 
ticular choice of cluster models in this paper follows the pa- 
rameters adopted by Chernoff fc Weinberg (1990) (see also 
Sigurdsson & Phinney 1995). The initial particle masses, 
positions and velocities are taken from an explicit, maximal 
A'' r ealisation of an isotropic Michie-King model of a clus- 
ter 



King 1962 



Da Costa fc Freeman 1976 



Gunn & Griffin 



197!;), as described in §2.2. The cluster's evolution along an 
orbit in a three component rigid model of the Galaxy (de- 
scribed in §2.3) is followed using a self-consistent field (SCF) 
code to calculate the mutual interactions of stars in the clus- 
ter (see §2.4.1). The SCF approach cannot follow evolution 
due to close encounters of individual stars in the system, so 
we repeat one case with the effects of two-body encounters 
included via diffusion coefficients to illustrate their impor- 
tance (see §2.4.2). 

This study is intended to isolate the response of a cluster 
to tidal effects with the dual purposes of identifying observ- 
able signatures of tidal interactions and providing a set of 
control models to compare to future studies which will in- 
clude more of the physics that might influence the cluster's 
evolution. In particular, we make the following simplifica- 
tions: 

(i) The distribution of particle masses is chosen to rep- 
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resent an evolved stellar population with a turnoff at ~ 
0.8Mq. We assume that continued stellar evolution is slow 
compared to the dynamical timescales involved. This as- 
sumption will be relaxed in subsequent papers, with explicit 
stellar evolution done in tandem with the dynamical evolu- 
tion. 

(ii) The initial cluster models are truncated at some finite 
radius, rr, corresponding to the 'tidal radius' of the King 
model. However, vt is not set to the limiting radius expected 
along each orbit, as a key purpose of the study is to inves- 
tigate observable effects of tidal shocks and stripping of a 
cluster that is not in exact equilibrium with the tidal field. 
The chosen models and orbits cover a range of interaction 
strengths - from those that only survive a few orbits to those 
that could last for the lifetime of the Galaxy. The scenario 
implicit in these initial conditions is that a cluster is formed 
very compact, with most of its mass well inside its actual 
tidal radius. Mass loss during stellar evolution subsequent 
to the cluster's formation leads to expansion of the cluster, 
until the tidal radius is reached. Here we assume that there 
is negligible loss of stars due to tidal effects until late in the 
cluster's evolution, and that we can start with a relaxed, 
mass segregated cluster with an exact King model profile, 
and follow subsequent mass loss due to tidal effects. Clearly 
this is an approximation to the real physics, but in order to 
isolate the tidal effects we are exploring, the initial condi- 
tions of the cluster must be close to equilibrium or internal 
dynamical evolution will instead dominate. 

(iii) As noted above, the effects of two body encounters 
on the internal dynamical evolution have not been included 
in the majority of our models. 

In this paper we restrict detailed discussions to effects that 
are insensitive to these simplifications. For example, in §4 we 
use our models to test how a cluster's extra-tidal population 
can be used to determine its mass loss rate - we expect 
this to be valid because although the mass loss rate itself 
will depend on the internal dynamics of the system (which 
we do not model exactly), the characteristics of the extra- 
tidal population are determined only by the tidal potential 
(which we do model exactly) . Conversely, we do not attempt 
to arrive at quantitative conclusions in our discussion of the 
evolution of clusters along different orbits (§3) since these 
would be influenced by all three of the above simplifications. 

2.2 Initial cluster models 

Four representative cluster models (labelled 0-3) are consid- 
ered, with properties summarised in Table |^. Model Oa was 
repeated with the same number density profile, but larger 
mass (Model Ob), and with different mass function indices 
(Models Oc and Od). Figure |l| shows the number density and 
mass density as a function of spherical radius, and surface 
brightness as a function of projected radius for each initial 
distribution. Masses are converted to visual magnitudes and 
luminosities in this and subsequent plots using Bergbusch & 
VandenBerg's (1992) isochrone for a 14 Gyr, [Fe/H]=-1.66 
cluster. Figure ^ illustrates the level of mass-segregation 
by showing the cumulative fraction of stars within a given 
radius for each mass group in the models. 

The cluster is generated as a Monte Carlo realisation of 
the Michie-King distribution function, using 8 discrete mass 



Figure 1. Initial number density (top panels), mass density 
(middle panels), and surface density (bottom panels) for each of 
the initial models. 
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Figure 2. Cumulative distribution of stars for each mass group 
in the labelled models. The highest/lowest curves correspond to 
the largest /smallest mass group (i.e. bin number ct increasing 
upwards). The mass range encompassed by a curve representing 
mass bin a can be found by multiplying the mass of the model 
(column 2 of Table |^ by the mass fraction dvon in bin a for the 
appropriate mass function index x (see Table bh . 
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Table 1. Description of Models. 



Model 


mass 


pc 


N 


Wq 


# density 


X 


'^dyn 

lO'^ years 


10® years 


orbits 


Oa 


0.34284 


10.6 


142365 


4 


50 


1.35 


6.19 


2.68 


pl-3,dl-2 


Ob 


2.73969 


21.2 


1136640 


4 


50 


1.35 


6.19 


17.66 


pl 


Oc 


0.20614 


6.24 


29904 


4 


50 


0.00 


3.61 


0.383 


Pl 


Od 


0.43807 


13.4 


269366 


4 


50 


2.50 


7.74 


6.00 


pl 


1 


1.22626 


7.36 


508669 


6 


1000 


1.35 


1.89 


2.62 


pl,p2,dl 


2 


2.86636 


14.1 


1186585 


9 


1000 


1.35 


3.28 


9.92 


pl 


3 


3.30330 


4.62 


1368230 


12 


1.0 xlO^ 


1.35 


0.573 


1.98 


pl,p2 



Columns: (1) Model number; (2) mass; (3) half mass radius; (4) number of stars; (5) King model; (6) central number density; 
(7) mass function index; (8) internal dynamical time (see equation |b[); (9) half-mass relaxation time; (10) orbits simulated. 



groups and a truncated, evolved power law initial mass func- 
tion as described in Sigurdsson & Phinney (1995). The struc- 
ture of each cluster is determined from a set of input parame- 
ters: the depth of the central potential, Wo, the central num- 
ber density, no, and the mean central dispersion, cr. Some of 
the observable structure parameters, such as the concentra- 
tion, core radius and hght profile, are derived quantities and 
depend on the mass function and mass-luminosity relation 
chosen. 

The particles representing the stars in mass group a(= 
1 ... 8), of mass m^, are assumed to have some initial distri- 
bution function, (7a(e), where e = — $ — is the energy 
and V is the velocity of the particle (both in the cluster 
centre-of-mass frame), and $ is the cluster potential. The 
distribution function is then given by 



(27rac.)3/2 



(3) 



where Gc is the core dispersion of mass group a and the no„ 
are the normalised densities of each mass group. The no„ 
are solved iteratively given the cluster structure parameters, 
according to the scheme described in Sigurdsson & Phinney 
(1995). 

For each mass function index x considered in this paper. 
Table ^ gives the average mass of stars and total mass 
fractior 



I'lUa assiguyd Lu bin u, along with the fiacLiun Ja of 
stars in' each bin that arc luminous. Stars that have evolved 
beyond the main-sequence turnoff are assigned to remnant 
stellar classes. In particular, stars with initial masses be- 
tween the turnoff and A.IMq axe assumed to have evolved 
to white dwarfs, and stars with initial masses between 8 
and 15M0 are assumed to have formed IAMq neutron stars. 
This follows the prescription in Chernoff & Weinberg (1990). 
Future work will consider other prescriptions for assigning 
remnant masses as a function of initial stellar mass. The 
presence of a substantial population of intermediate mass 
white dwarfs changes the present day mass function from 
the simple power law of the initial mass function. In partic- 
ular, there is a pronounced 'bump', an excess of stars with 
masses ~ Q.^Mq in the distribution function, which compli- 
cates interpretation of the local luminosity function in terms 
of the mass function, for intermediate mass stars. 

Each cluster is generated by a simple acceptance- 
rejection algorithm, selecting particle masses, positions and 
velocities from the distribution function. The masses are as- 
signed in discrete mass bins, with the value of the mass in 



each bin being the mean mass of that interval, weighted by 
the evolved mass function. That is, the mean mass in each 
bin allows for the presence of stellar remnants in that bin. 
Future models may incorporate a continuous distribution of 
stellar masses. The positions are drawn from a radial grid 
(typically of about 200 points, spaced to sample the density 
profile efficiently). A particle selected to be at radius is 
assigned to some radius ± e where e is a uniform random 
variable on the grid space interval. The initial distribution 
is thus a set of thin stepped radial shells, with local density 
that deviates slightly from the true density profile. Phase 
mixing ensures that the model settles down to a smooth 
representation of the cluster on a dynamical timescale. The 
initial model has a virial ratio ^ 1 ± 1/ VN , and is station- 
ary and stable. Cluster models were run in isolation using 
the SCF code (see §2.4.1 and Hernquist & Ostriker 1992): 
the initial mass segregation is robust and there is no spu- 
rious evolution. The mass segregation corresponds to ther- 
mal equilibrium for the cluster structure parameters (cen- 
tral potential and initial mass function) chosen. This is the 
equilibrium state expected for the internal distribution of 
masses within the cluster once the time scale for mass loss 
due to stellar evolution becomes longer than the core relax- 
ation time, and provided that tidal effects have not yet had 
a strong effect on the cluster. There is also some e vidence 
for primord i al mass segregati on in young clusters (Fischer 
et al. 1998; Elson et al. 1998), possibly due to bias in the 



formation of massive stars towards high density regions, or 
due to rapid relaxation in the young cluster. Mass segrega- 
tion continues on a relaxation time scale, and would lead to 
core collapse in the absence of other effects. 



2.3 Cluster orbits in the Galactic potential 

A three-component model is used to represent the Galactic 
potential: 'I' = 'I'disk + 'I'sphcr + ^'haio, in whi ch the disk is 



represented by a Miyamoto-Nagai potential (Miyamoto & 



Nagai 1975 ), the spheroid by a Hernquist potential (Hern- 
quist 199C ) , and the halo by a logarithmic potential: 



GAfdis 



sphcr 



+ (a + V-z^ + h'^Y 
r + c 



*haio = fhalo ln(r -)- d ). 



(4) 

(5) 
(6) 
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Table 2. Adopted Mass Functions with Index x 



a 




X = 0.0 






X = 1.35 






X = 2.50 








drria 




rria 


dma 






drria 


U 


1 


0.12629 


0.01847 


1.00000 


0.12346 


0.23393 


1.00000 


0.12114 


0.50252 


1.00000 


2 


0.22461 


0.04963 


1.00000 


0.21330 


0.29045 


1.00000 


0.20449 


0.33330 


1.00000 


3 


0.34806 


0.02601 


1.00000 


0.34600 


0.08360 


1.00000 


0.34427 


0.05484 


1.00000 


4 


0.44265 


0.03596 


1.00000 


0.43956 


0.08355 


1.00000 


0.43695 


0.04163 


1.00000 


5 


0.57392 


0.12309 


0.34155 


0.56679 


0.13430 


0.52623 


0.56126 


0.03946 


0.67655 


6 


0.71009 


0.16844 


0.32948 


0.70416 


0.11552 


0.58798 


0.70078 


0.02525 


0.77595 


7 


0.99995 


0.22745 


0.00000 


0.96585 


0.04281 


0.00000 


0.93946 


0.00264 


0.00000 


8 


1.38462 


0.35095 


0.00000 


1.36336 


0.01583 


0.00000 


1.34231 


0.00036 


0.00000 



Columns: (1) — mass group; (2), (5) and (8) - average mass m^; (3), (6) and (9) - mass fraction dma; (4), (7) and (10) — 
fraction of luminous stars assigned to bin a. 
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Figure 3. Locus of polar orbits pi (left hand panel), p2 (middle 
panel) and p3 (right hand panel). The Galactic plane is at Z = 0. 
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Figure 4. Locus of disk orbits dl (left hand panel) and d2 (right 
hand panel). The solid line is for motion perpendicular to the 
Galactic plane (i.e. Z — X motion) and the dotted line is for 
motion in the Galactic plane (i.e. Y — X motion). 



Here, Mdisk = 1.0 x 10'\M,pher = 3.4 x 10'°,Uhaio = 
128, a — 6.5, b = 0.26, c = 0.7, and d — 12.0, where masses 
are in Mq, velocities are in km s~^ and lengths are in kpc. 
This choice of parameters provides a nearly flat rotation 
curve between 1 and 30 kpc and a disk scale height of 0.2 
kpc. The radial dependence of the z epicyclic frequency (k^) 
in the disk between radii at 3 and 20 kpc is similar to that 
of an exponential disk with a 4 kpc scale length. 

To contrast the evolution of the models in different tidal 
fields we choose orbits roughly 3kpc, 5kpc and lOkpc from 
the Galactic centre, three with polar orientations (pi, p2 



and p3) and two near the Galactic disk (dl and d2). Figures 
^ and ^ show the paths of these orbits and Table ^ gives 
their peri- and apo-galacticon and radial time periods. 

Hereafter we refer to the simulation representing the 
evolution of cluster Model m along orbit o as Model (m, o). 



2.4 Integration methods 

Individual particle trajectories are integrated using a 
leapfrog integration scheme, with accelerations calculated 
from 



r = V(*(r) + <I>(r)) + ajyf + aMck, 



(7) 



where the Milky Way is represented by the rigid potential 
given in §2.2, and the cluster's internal potential <1> is 
calculated from the particles' masses and positions using the 
SCF code described in §2.4.1. The effect of accelerations due 
to two-body effects (i.e. terms adyf -I- akick) were included 
in only one simulation using the method described in §2.4.2. 
Dynamical friction of the cluster's orbit was ignored. In all 
cases, the time step was smoothly varied along the orbit 
between Tdisk/lOO and Tdyn/lOO to ensure that energy was 
conserved to better than 1 percent of the initial internal 
potential energy of the cluster. Here 



(8) 



is the timescale for a disk passage, where b is the vertical 
scale of the Galactic disk (see equation [Q), and Vz is the 
z-component of the satellite's velocity as it crosses the disk 
plane at Z — 0, and 



Idyn 



2Gme, 



(9) 



is the internal dynamical timesca le for a system of total mas s 
JTisat and half- mass radius rhait (Binney & Tremaine 1987). 
The timescales for the orbits and models are given in Tables 
^ and 1^ respectively. 



2.4-1 Self-consistent field code 

The SCF code uses a bi-orthogonal basis function expan- 
sion to calculate the internal potential of the cluster from 
the individual particle positions and masses (Hcrnquist & 
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Table 3. Description of Orbits. 



Orbit 


^pcri 


^apo 




^disk 


Tdisk 


^bulgc 


-^bulgc 




kpc 


kpc 


10^ years 


(km s~-'-)/kpc/Mycar 


IQ^years 


(km s~-'^)/kpc/Myear 


lO^years 


Pl 


1.1 


3.5 


8.0 


30 


1.0 


20-30 


3. 


p2 


3.0 


5.8 


14.5 


20 


1.0 


2-3 


10. 


p3 


9.5 


12.5 


32.0 


5-10 


1.0 


0.5 


20. 


dl 


2.9 


3.15 


8.5 


30 


2.0 


3-4 


12. 


d2 


4.6 


5.4 


14.5 


30 


2.2 


1 


20. 



Columns: (1) orbit; (2) closest approach; (3) apocentre; (4) azimuthal orbital time period; (5) defined in equation (p^; (6) 
timescale for disk passage (see equation Isj); (7) and (8) as columns (5) and (6) but for bulge passages. 



Ostriker 1992). The expansion is sensitive to large-scale fluc- 
tuations in the cluster's potential but smooths local poten- 
tial fluctuations arising from the discrete particles. Hence 
the SCF approach will underestimate relaxation when (as 
in our case) the number of particles in the simulation is 
equivalent to the number of stars in the system represented. 
The computations using this scheme (i.e. not including two- 
body encounters) were performed at the National Center for 
Supercomputing Applications (NCSA) with a version of the 
SCF code which h ad been parallelised t o run on the Con- 
nection Machine 5 ( Hernquist et al. 1995 ). On average, each 
step required 4 x 10~ cpusec/particle/processor. 



2.4-2 Two-body relaxation calculation 

In order to explore the consequences of continuing mass seg- 
regation and two-body relaxation on the structure of the 
cluster, we incorporated a Fokker-Planck diffusion scheme 
i nto a version of the SCF code parallelised to run on the T3E 
(Sigurdsson et al. 1997) at the Pittsburgh Supercomputing 
Center. The results of Model (0a,p3), run with and without 
diffusion are compared in in §3.3. 

The Fokker-Planck scheme calculates the flrst and sec- 
ond order velocity diffusion coefficients, D{Avi), D{AviAvj) 
(Binney & Tremaine 1987). The diffusion is directly incor- 



porated into the explicit time evolution of the particles by 
including the last two terms in equation (^), where adyf is 
the dynamical friction experienced by the star, and akick 
is the effective acceleration due to scattering by individual 
stars in the cluster (see Sigurdsson & Phinney 1995 for dis- 
cussion) . 

To calculate adyf and akick, we choose an orthonor- 
mal basis local to each particle defined by the particle's ve- 
locity V and position x in the cluster frame. This gives us 
three independent diffusion coefficients, I?(Aii||), D(Aiiy), 
and D{Av'^). To calculate those we need the local d ensity. 



which is obtained dir ectly from the SCF expansion (Hern 
quist & Ostriker 1992), and a local velocity distribution. We 
approximate the velocity distribution as a Gaussian, and 
calculate it every N^, ~ 100 integration steps, by sampling 
the velocity distribution in radial bins, calculating the lo- 
cal dispersion, and using spline interpolation to obtain the 
dispersion as a function of radius. Provided the dynamical 
evolution of the cluster is slow compared to the dynami- 
cal time scale, the sparse updating of the dispersion profile 
is adequate, and provided the cluster is close to being re- 
laxed, the approximation of taking the velocity distribution 
as an isotropic Gaussian is acceptable. The diffusion ap- 



proximation incorporates a Coulomb logarithm term whose 
magnitude is uncertain to a magnitude comparable to other 
errors incurred by the approximations we make. Near the 
truncation radius, the interpolation of the dispersion must 
be positive definite or negative dispersion may be inferred. 
Outside the truncation radius the diffusion coefficients may 
be set to zero, or the local dispersion approximated by the 
Keplerian velocity due to the total enclosed mass. 

Given the diffusion coefficients, then adyf — D{Av^^) 
and we model akick by random fiuctuations in velocity, Av, 
where akick = -§7 , with 

Avl =if{D{Avf^)At) (10) 
,i{DiAvi)At), (11) 



Av'i 



where q is a random number with zero mean and unit stan- 
dard deviation, chosen here from a normal distribution. 

This choice of diffusion coefficients provides a fast and 
surprisingly good approximation to two-body relaxation in 
clusters, and allows a fair representation of the effects of 
two-body relaxation on the structure of the cluster in the 
presence of other perturbing dynamical processes. For iso- 
lated clusters, this scheme provides a realisation of evolution 
towards core collapse on relaxation time scales, and can fol- 
low the evolution of the cluster over several orders of mag- 
nitude in central density. Ultimately though, with a finite 
number of expansion terms, the SCF code fails to resolve 
the resulting density cusp and core collapse to infinite den- 
sity is not observed. Details of the scheme will be discussed 
in another paper. 

2.5 Analysis methods 

2. 5. 1 Isolating the bound population of the cluster 

In parts of this study, identical analyses are performed on 
the particle data from the simulations both with and without 
the unbound particles, and the results are compared. Much 
of the discussion of the interpretation of observations rests 
upon the identification of stars that are cluster members. We 
define members of the cluster to be the maximum set of mu- 
tually bound particles. This set is found iteratively, starting 
from the set of all particles, by: (i) calculating the internal 
potential field of the particles in the set considered; (ii) find- 
ing the velocity of each particle with respect to the velocity 
of the minimum of this potential field; (iii) defining each 
particle's internal kinetic energy to be the kinetic energy of 
this relative motion; (iv) labeling those particles whose in- 
ternal kinetic energy is greater than their internal potential 
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energy as 'unbound', removing them from the set considered 
and returning to step (i). Steps (i)-(iv) are repeated until no 
new particles in an iteration are labelled 'unbound'. Using 
this method we find that the instantaneous mass bound to 
the cluster is a monotonically decreasing function of time, 
and the number of stars that later become bound again to 
the cluster once they have first been classified as unbound 
is negligible. 

2.5.2 Extrapolating mass functions 

Each particle in our simulations is assigned the average mass 
nia of the mass bin a that it occupies. The first four mass 
bins contain only luminous matter, the next two contain a 
fraction fa of luminous stars and the two most massive bins 
contain only stellar remnants (see Table ^). Hence, we can 
trivially find the mass function of observable stars in any 
sample at these average points by calculating 

fdN\ _ fgNg 
\dmj a Ania ' 

where Arria is the width of the mass bin and Nc, is the 
number of stars in the bin. Since mass functions are typically 
calculated from stars near the turnoff in a globular cluster, 
we estimate x using only the two heaviest luminous mass 
bins to be 



(12) 



LOGjNeffiAms) - LOGiNsfsAme) 
LOG{me) - LOGims) 



^-1 



-3 



□ Orbit p 1 




1 1 1,, 1 1 


A Orbit p2 
Q Orbit p3 


(Oa:pl)^ 


□ (papi) 


a Orbit dl 
i Orbit d2 


yn(a:pl) 


,(Oa:dl) 




A(Oa:pS) 


A(0a:d2) 


□■(l:pl) 

□ (3:pl) 






«(3:dl) 




- 




.(l:dl) 




A(l:p3) 

O 

(0a:p3) 
< 1 < < , , 




, , 1 < , 


-3 


-2 


-1 



LOG [(AE,„^ >,,,,/a^] 

Figure 5. Average potential change per disk shock for each model 
as a function of shock strength, in units of the clusters' internal 
velocity dispersion a. Each point is labelled with the model and 
orbit. The dotted line indicates the trend expected when the im- 
pulse approximation is valid. 



(13) 



3 GENERAL EVOLUTION IN A TIDAL FIELD 
3.1 Cluster heating 

We expect mass loss to occur predominantly where the tidal 
field of the Milky Way is strongest - i.e. at pericentric points 
along the orbit and during disk passages. In column 5 of 
Table |^ we give values for 



(14) 



For a globular cluster whose physical scale is given by its 
half mass radius rhaif, this can be used to estimate the spe- 
cific energy change due to the disk passage in the impulsive 
regime, 

1 



AEi„ 



(^diskrdisk'"ha 



(15) 



(Bpitzer 1958) which we shall subsequently refer to as the 
shock strength, where Tdisk is the timescale for the disk pas- 
sage (see equation |^). We expect this estimate to be ap- 
propriate when Tdyn 3> Tdisk, and can otherwise use it as an 
upper limit on the average energy change. Columns 7 and 8 
in Table ^ give similar estimates for the importance of bulge 
shocking for these orbits (with Abuige ~ d^*&buigc/rf^^), 



which ih cleaily LuuipeLili 



ith disk shuLkiug fui uibii pi 



(see uilai. Hut &l OiLiikti 1000 fui a dijtujjiuu uf thij 
effect). ' In our choocn Galactio potential, all thcac orbits lie 
within the core radius d = 12kpc of the halo component (see 
equation |Q), so its tidal field is roughly constant along the 
orbit and halo shocking is unimportant (in contrast, this 
is likely to be an important factor in the evolution of the 
Sagittarius dwarf spheroidal galaxy - see Johnston, Spergel 
& Hernquist 1995). 



To test the accuracy of the impulse approximation, we 
calculated the change in the cluster's internal potential en- 
ergy per unit mass (A<1>) and the impulsive estimate for this 
change (AiJimp - from equations and |^) for each 

of the A'^disk disk passages in each simulation. The points 
in Figure ^ show the average values of these quantities in 
the simulations, in units of the cluster's velocity dispersion, 
a^. The models in the bottom left hand corner of the plot 
are expected to be the most resilient to the Milky Way's 
infiuence. 

There is a large spread in the response of the clusters' 
potential energy for a given shock strength, AiSimp. The 
clusters on disk orbits (filled symbols) are shocked less no- 
ticeably than those on polar orbits (open symbols), because 
the timescales for the disk passages in the former case are 
longer, and the impulse approximation will only be valid in 
the outer region of the cluster. In general. Figure ^ shows 
that while this crude estimate serves as a useful guide to 
the importance of shock heating for a cluster in a given or- 
bit a much more detailed calculation of the interaction of 
the internal and external dynamics is needed for an accu- 
rate assessment of cluster evolution in general. More rigorous 
analytic estimates for tidal heating have been discussed ex- 
tensively elsewhe r e in the literature ( Weinberg 1994 ^ Kundic 
Gnedin & Ostriker 1997; pohnston ct al 



Ostriker 1995 
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3.2 Differential mass loss 

Figures I and@ follow the bound mass fraction remaining 
as a function of time for the first 0.5 Gyrs of evolution along 
the polar and disk orbits. On each of these figures we show 
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time (Gyears) 

Figure 6. Bound mass fraction as a function of time for the mod- 
els along the polar orbits pi (top panel), p2 (middle panel) and 
p3 (lower panel). The dashed lines label the time of disk passages 
and the dotted lines label the time of pericentric passages. 
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Figure 8. Bound mass fraction, limiting and half mass radii for 
each mass group as a function of time for Model (l,pl) (left hand 
panels) and Model (0a,p3). The vertical dashed and dotted lines 
give the time of disk and pericentric passages respectively. The 
outer solid lines in each panel are for the least and most massive 
stars and the dotted lines are for the intermediate mass stars. 
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Figure 7. Bound mass fraction as a function of time for the mod- 
els along the disk orbits dl (top panel) and d2 (bottom panel). 
The dashed lines label the time of disk passages. 



the time of disk passages as vertical dashed lines, and the 
pericentric points along the orbits as vertical dotted lines. As 
noted in the previous section, the mass loss rate is greatest 
at these locations along the orbit. A comparison of these 
plots with Figure |^ confirms the loose correlation between 
mass loss rate and shock strength. 

Figure ^ compares the effect of the tidal field on each 
individual mass group in Model (l,pl) (left hand panels) 
with those in Model (0a,p3) (right hand panels). The solid 
lines are for the smallest and largest masses and the dotted 
lines are for the intermediate ones. The top panels illustrate 
differential mass loss due to mass segregation, as the most 



weakly bound (and therefore lower mass) stars are preferen- 
tially stripped. The lower two panels give a measure of the 
response in the spatial distributions - is calculated as the 
average radius of the outermost 1 percent of bound particles 
in each group, and rhaif is calculated as the radius contain- 
ing half the mass of each group. Again, the compressive disk 
shocks (and subsequent expansion due to the energy input) 
are clearly seen. The effect is more dramatic in the left hand 
panels since this model is more susceptible to the influence 
of tides (as demonstrated by its position in Figure |^, and 
large mass loss rate). Similar characteristics are seen in all 
the models. 

Figure ^ repeats Figure^ for Models (Oa,pl) (left hand 
panels) and (Ob, pi) (right hand panels). These models had 
identical mass density profiles, but with different radial 
scales (and hence total masses) . The comparison of the two 
demonstrates that the general characteristics of the response 
for models with the same density profiles subjected to the 
same tidal field are identical. 

Figures ^ - ^ illustrate our intuitive understanding of 
tidal mass loss in a mass-segregated system. As an example, 
the solid lines in Figure hd show the bound mass fraction 
for each mass group as a function of time for Model (Oa,dl). 
The dashed lines show the prediction if the initial model 
were simply truncated in radius at the surface enclosing the 
same mass as is instantaneously bound to the system in the 
simulation. The dotted lines show the prediction if the initial 
model is instead truncated in binding energy. As might be 
anticipated, initial binding energy provides the more accu- 
rate of these simple guides to which stars are lost first from 
the system. This behaviour is seen in all our simulations, as 
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Figure 9. As Figure H but for Models (Oa,pl) and (Ob, pi). 
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Figure 10. The solid lines show the evolution of the bound mass 
fraction for each mass group in Model (Oa,dl). The dashed lines 
show the fraction of each mass group within a spherical surface 
in the initial model which contains the same total mass as is 
instantaneously bound to the satellite. The dotted lines show the 
fraction of each mass group within an energy surface in each initial 
model which contains the same total mass as is instantaneously 
bound to the satellite. 



Figure 11. The solid lines show the evolution of the mass func- 
tion index x along each of the polar orbits. The dashed lines give 
X for the mass within a spherical surface in each initial model 
equivalent to the bound mass remaining in the simulations at 
that time. The dotted lines show x within an equivalent energy 
surface in each initial model. 
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Figure 12. As Figure |ll| but for the disk orbits. 



shown in Figures |l^ and The solid lines show the global 
mass function index function of time for each simu- 

lation, the dashed lines show the result of estimating x by 
truncating the initial model in radius and the dotted lines 
show X for truncation in energy. Note in general that the 
agreement is good, except when there has been substantial 
evolution of the system (e.g Models (Oa,pl), (Oc,pl)). 

In summary, the results in this section demonstrate 
that: (i) the impulse approximation provides a rough guide 
to the mass loss rate from a system; (ii) tidal mass loss 
along a given orbit is determined by the density profile of 
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Figure 13. Bound mass fraction as a function of time for the 
labelled mass groups in Model (0a,p3) run with (dotted lines) 
and without (solid lines) diffusion. 
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Figure 14. Comparison of initial (stars) and final properties 
of Model (0a,p3) run with (open squares) and without (closed 
squares) diffusion. 



the cluster; and (iii) tidal shocking of a mass-segregated sys- 
tem leads to evolution of the mass function, which can ap- 
proximately be predicted for a given mass lost by simply 
truncat ing the initial distribution in energv. Clearlv these 



general [trends fit in with theoretical interpretations of ob- 
s erved correlations of c luster properties with Rgc and Z 
(gapaccioli et al. 1991; gapaccioli, Piotto fc Stiavelli 1993 



Djorgovski^'^iotto'TrCapa^^^ However, a quantita- 

tive comparison of our simulations with observations would 
require more detailed analytic models both of the tidal ef- 
fects seen in the simulations and of other physical processes 
not included (one example is given in the next section) and 
is beyond the scope of the current study. 



3.3 Competition of tides and relaxation 

Figure ^ compares the mass loss rate for each group in 
Model (0a,p3) evolving over 8 Gyrs both with (dotted lines) 
and without (solid lines) diffusion included. The half-mass 
relaxation time for this cluster is Treiax ~ 2.68 Gyrs (see 
Table [l|). The evaporatio n time for a cluster is ~ 100 Trdax 
(Meylan & Heggie 1997), so we expect only a few percent 



change in bound mass fraction due to relaxation over the 
course of the simulation, as is seen in the figure. Relaxation 
contributes to mass loss by increasing the number of low 
mass stars lost and decreasing the number of high mass stars 
lost, with a net increase in mass loss. Hence, the evolution 
of the mass function will be modified from the purely tidal 
response described in the previous section. Figure [Til shows 
that the density and velocity properties of the model evolved 
with (solid squares) and without relaxation (open squares) 
are similar (the initial model is shown in stars). 

These figures demonstrate that our neglect of relaxation 
in the models considered in this paper will not substantially 



affect their evolution since a significant fraction of mass is 
lost due to tidal infiuences over a few relaxation times. In 
contrast, Pal 5 is an example of a cluster whose evolution 



is th ought to be dominated by relaxation (Rosenberg et al 



199S ) and a natural extension to our work would be to con- 
sider models whe re the tidal and evaporati on mass loss rates 
are comparable (Vesperini & Heggie 1997). 



4 INTERPRETING OBSERVATIONS OF 
TIDALLY DISRUPTING, 
MASS-SEGREGATED SYSTEMS 

4.1 'Observing' the simulations 

Observed properties of globular clusters are often inter- 
preted as resulting from the tidal influence of the Milky 
Way. In this section, we 'observe' our models to confirm the 
validity of some common interpretations and explore other 
signatures of tidal disruption that future studies might be 
sensitive to. We focus our analysis on Model (0a,p3) as an 
example of a cluster with a steady mass loss rate that could 
survive for the lifetime of the Galaxy. We examine the ob- 
servable properties of this model at the end of the simulation 
(after 8 Gyrs of evolution) , and compare and contrast it with 
the same model at earlier times and with the other models 
in our sample. 

Since the Sun's position along the Solar Circle in the 
simulations is arbitrary, rather than assuming a single view- 
point we examine each set of particle data along three per- 
pendicular axes defined by the cluster's instantaneous po- 
sition and velocity: the x axis lies along the vector joining 
the cluster to the Galactic centre; the z axis is perpendicu- 
lar to the cluster's position and velocity vector (i.e. looking 
down on the orbital plane); and the y axis is perpendicular 
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Figure 15. Number count surface density profiles along labelled 
views (see §4.1) and down to various limiting magnitudes {Mto 
is the magnitude of the turnoff) for the final bound population 
(filled squares) and all stars (open squares) for Model (0a,p3). 
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Figure 16. Number count surface density profiles along the x- 
axis for the bound population (filled squares) and all stars (open 
squares) at random points along the orbit for labelled models. 
The dashed lines show the surface density for extra-tidal material 
from equation ( p^ ) and the dotted lines show the adopted rtrcalj 
in each case, for use in the mass loss estimates. 



to these two, along the orbit of the satellite (i.e. it coincides 
with the velocity vector when the Galactocentric radial ve- 
locity is zero) . We also make the simplification that the lines 
of sight across the face of the cluster are perfectly parallel 
(i.e. the observer is at infinite distance from the cluster) and 
state distances in our projected coordinates in pc rather than 
arcmin. This will not significantly affect our analysis of den- 
sity profiles in §4.2, but can alter the measurement of cluster 
rotation, and we discuss this simplification further in §4.3. 

4.2 Star count profiles 

4-2.1 Qualitative interpretation 

Figure |l5| shows the number surface density for all stars (top 
panels), stars brighter than 4 magnitudes below the turnoff 
(middle panels) and stars down to the turnoff (bottom pan- 
els) at the end of the simulation of Model (0a,p3) with the 
three labelled viewpoints along the axes defined in §4.1. The 
closed symbols are for stars still bound to the satellite and 
the open symbols include unbound stars in the calculation 
(where the 'bound' and 'unbound' populations are separated 
following the method described in §2.5.1). 

In this Figure, there is a clear break in the slope of the 
profile defined by the open symbols, approximately corre- 
sponding to the point where the closed and open symbols 
become distinct and the fraction of unbound stars in an 
annulus becomes significant. Such 'extra-tidal' stars have 
been detected in several studies of dwarf spheroidal galaxies 



(|IH95|; [Kuhn et al. 1996D and globular clusters (|G95| prill 
mair 1998 ). Our investigation confirms the common inter- 
pretation that these can be identified as stars escaping from 
the satellites. Figure |l^ demonstrates that this is a general 
result by repeating the x-axis view for all stars in Model 



(0a,3p) at earlier times (left hand panels) and at random 
points for several of the other models (as labelled). 

4-2.2 Quantitative interpretation 

Tremaine (1993) pointed out that the change in the or- 
bital frequency of a star torn from the satellite at ra- 
dius rbreak (the poiut at which the slope of the surface 
density profile changes) should approximately be given by 
Af2 ~ rhieakdQ/dR, where fl{R) is the frequency of a cir- 
cular orbit in the parent galaxy at radius R- Equivalently, 
the orbital energies in tidal debris from a satellite orbiting 
in a potential <1> will be spread over a characteristic range 
e = rtidcd^ / dR. Johnston (1998) tested this simple physical 
argument by looking at the spread in energies in streamers 
seen in numerical simulations of tidal disruption and found 
that e < \AE\ < 2e. Hence, debris will spread over an angu- 
lar distance comparable to the size of the cluster (rbroak/^) 
in a time (rbroak/i?)/2Ar2 ~ Torb/iT, where Torb is the az- 
imuthal time period of the orbit. This suggests that we can 
estimate the average surface density of stars in an annulus 
between rbreak and r from the centre of the cluster to be 



dm (r — rbroak) To^h 



dt T^brcak 
2 



(16) 



/[n{r 



^brcak / 



where dm/dt is the mass loss rate from the cluster. The func- 
tion g depends on the angle 9 of our line of sight with the 
plane perpendicular to the direction of motion of the satel- 
lite. Eventually, the debris spreads out along the satellite's 
orbit so we can take g{9) — cosO as a crude approximation. 
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Figure 17. Known mass loss rate df /dt (from Figures Q and ^ 
for the labelled models plotted against the estimates (df /dt)i (left 
hand panel) and {df/dt)2 (right hand panel). The solid squares, 
open squares and stars show the results when viewed from along 
the X, y and z-axis respectively. The open triangles repeat the 
y-axis view but with g{9) = 1. The solid line shows where the 
estimated and known mass loss rates agree, and the dotted lines 
are for factor of 2 discrepancies. 



In fact, mass tends to leave the satellite at the inner and 
outer Lagrange points along the vector joining the satellite 
to the centre of the galaxy, so the geometry of the streamers 
close to the satellite is rather more complicated than this. 

We can differentiate equation ( [l6| ) to find an approxi- 
mate expression for the absolute surface density 

1 dm Torb 1 



(17) 



g{e) dt TT 27rrbroakr' 

This estimate is overlaid in dashed lines on the profiles 
shown in Figure |l6| using dra/dt averaged over each simula- 
tion, Torb from Table ^ and rbrcak indicated by the vertical 
dotted lines in each panel. The open squares and dashed 
line agree well in all cases except for Model (3,pl). In this 
model, the simulation was run for only a few orbital periods 
so the debris only had time to disperse a few rbrcak irom the 
satellite and the streamers were not well populated. 

Equation ( p^ ) and Figure |l6| demonstrate that we ex- 
pect the extra-tidal population around a cluster to have 
Sxt(?') cx . In observations of real globulars, the sur- 
face density of extra-tidal stars has been found to fall as 
Sxt(r) ~ r"^ with -5 < 7 < -0.7 (G95; Zaggia, Piotto & 
Capaccioli 1997). This suggests that the measured 7 might 
be used as a test of whether the globular is sufficiently ob- 
scured, either by tidal debris further along its own tidal tail 
or by the Galactic field, that the above interpretation is in- 
valid. 

Note that these estimates are independent of the 
method of mass loss from the cluster since the dynamics 
once the star is lost is determined only by the external field, 
and hence should apply equally to mass lost by tidal strip- 
ping, shocking, or evaporation due to relaxation. 

4-2.3 Estimating mass loss rates using extra-tidal stars 

Equations ( |l^ ) and ( [l7|) suggest two ways of estimating the 
fractional mass loss rate [df /dt) from a cluster using obser- 
vations of extra-tidal stars. 

If the extra-tidal population is well defined out to a 



radius r-xt with Ext ~ we can count the number of stars 
"-break within rbrcak (the point where there is a break in 
the slope of the surface density profile) and the number of 
extra-tidal stars nxt between rbreak and rxt and use them to 
find 



^brcak 



nxt 



^xt ^brcak ^break '^orb 



(18) 



Since the timescale for stars to diffuse beyond a few rbrcak 
is an orbital timescale, so long as rxt > 2rbreak this estimate 
should be sensitive to the average mass loss rate rather than 
the instantaneous one. The left hand panel of Figure [l^ plots 
the known fractional mass loss rate, df /dt, for each of the 
models with fractional mass loss rates less than lGyr~^ (see 
Figures ^ and |^ against [df /dt)\ calculated from equation 
([l8|). To find {df /dt)i, we evaluated the ratio nxt/wbroak for 
rxt = 2rbrcak, rxt = 5rbroak and rxt = lOrbrcak- This pro- 
cess was repeated looking along the a::-axis (solid squares), 
2/-axis (open squares) and z-axis (stars). The open triangles 
show the estimate made along y-axis if the viewing angle 
is not taken into account (i.e. g{9) = 1). The solid line in 
the figure shows where the estimated and known mass loss 
rates agree, and the dotted lines are for factor of 2 discrep- 
ancies. The figure suggests that we can use this technique 
to estimate the mass loss rate from a Galactic satellite pro- 
vided our viewpoint lies close to the plane perpendicular to 
the satellite's velocity vector (equivalent to the xz-pla,ne in 
our projection). Otherwise, we are likely to poorly estimate 
the rate of mass loss as the debris along the satellite's orbit 
confuses the calculation. 

If the extra-tidal population is either not well defined 
or non-existent we can estimate {df/dt)2 as an upper limit 
for the mass loss rate from Ext (rbreak), where rbreak is taken 
to be the point either where the slope of the profile changes 
or the last point where the surface density is separable from 
the background. Then, from equation (|l7[), 



\dt / 2 "break 



1" r, 2 

ZTrrbrcak? 



(19) 



which is shown in the right hand panel of Figure |l^. As with 
the first case, the view perpendicular to the velocity vector 
provides the best estimate. 



4.3 Line of sight velocity distributions 

Figure |l^ summarises the line-of-sight velocity analysis at 
the end of the simulation for Model (0a,p3) from the stated 
viewpoints and with the observer assumed to be at infinite 
distance from the cluster (i.e with the individual lines of 
sight across the cluster assumed to be parallel to the one 
to the cluster centre - the effect of this simplification is 
discussed in §4.3.1). For each viewpoint the average Vios and 
dispersion aios in the velocities of stars are calculated in bins 
lying in strips along the two perpendicular axes across the 
cluster. Figure |l^ demonstrates the general nature of these 
plots, repeating the x-axis view of Figure ^ at different 
orbital phases of this model and for several different models. 
The solid symbols show the analysis for the bound stars, and 
the open symbols shows the results for all stars. 
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Figure 18. Line-of-sight velocity vi^^ and dispersion trios for 
bound (filled squares) and all stars (open squares) in Model 
(0a,p3). The line of sight for each column of panels is labelled 
in the top row (see §4.1 for axis definitions). 



■B, ° 



. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 IJ 

L {0,p3) i 


. 1 1 1 1 1 1 1 1 1 1 U 1 1 J 

^ (i,pi) °y ^ 

r 1 


. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 
L (l,p3) J 






□ ° 



y <p«) 



y (pc) 



^^.r^p^-^^SpI.^i^ 


- M 1 1 M 1 1 1 1 M U 

- □ □ n - 

A- ^ 

■ 1 1 1 1 1 T 1 1 1 1 1 1 1 1 1 - 


-1 M 1 1 1 M 1 1 1 1 M 1 - 

■ , ,°i *r , , 1 , , , ■ 1 , , ■ 


-100 100 -100 100 -100 100 

y (pc) y (pc) 7 (pc) 


. , , 1 1 1 1 1 1 1 , , , 1 , 1 _ 

(0,p3) 

- , , 1 , , , , 1 , , , , 1 , , = 


. 1 1 1 1 1 1 1 1 , 1 ^ , i^^.i 1 _ 
L (i.di) J 

^ / 1 

; , p° 1 , , , , 1 , , ; 


. . 1 1 1 1 . 1 1 . ^ . 1 1 1 1 . 

L. - ] 

- , , 1 , , , , 1 , , , , 1 , , - 


-100 □ 100 -100 100 -lao o loo 

y (pc) 7 (pc) y (pc) 


. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 _ 


°, ° 1 , , 1 , f , , 1 , , : 


^ s:= ^ 

□ ■ ■ □ 

- , R 1 , °, , 1 , . „ , 1 , , - 



-100 100 -100 100 -100 100 

y Cpc) y (pc) y (pc) 



Figure 19. Repetition of panels x-axis view of Figure |l^ for the 
labelled models at random points along the orbit. 



4-3.1 Average velocities and cluster rotation 

The average velocities are zero except for the viewpoints 
which are sensitive to the cluster's orbital motion. In the 
panels corresponding to the x-aixis view of Figures |l8| and 
|lS| , the dotted lines indicate the expected velocity gradient if 
the cluster were rigidly co-rotating with its orbit. The bound 



stars follow this line fairly closely, with some indication of 
a smaller gradient in velocities towards the centres of the 
clusters. The unbound stars appear to be rotating faster 
than this as they move to orbits with higher/lower angular 
velocities to form the leading/trailing streamers. 

In our analysis, two effects are clearly contributing to 
the measured rotation of the cluster - the intrinsic rotation 
of the cluster (in our simulations, roughly corresponding to 
co-rotation with the tidal field), and the velocity gradient 
in stripped material. Observers also have to contend with 
'perspective rotation' (Feast, Thackeray fc Wcssolink, 1961 



Merritt, Meylan fc Mayor 1997t [Prukicr ct al. 1997[ ) from the 
projection of the tangential velocity onto the non-parallel 
lines of sight across the cluster. If we were indeed viewing a 
rigidly co- rotating cluster from the centre of the Galaxy, the 
intrinsic and perspective rotation would cancel each other 
out and we would be sensitive only to tidal stripping, but 
if we were closer to/farther from the cluster than the cen- 
tre of the Galaxy, we would be dominated by the perspec- 
tive/intrinsic rotation. 

Unfortunately, we cannot assume that tidal torquing of 
a real satellite necessarily results in co-rotation. The rota- 
tion seen in our simulations is likely to be an artifact of the 
near-circular orbits we have chosen. In models, run along 
eccentric orbits, of the disruption of the Sagittarius dwarf 
galaxy, Johnston, Spergel & Hernquist (1995) found that 
the line of sight velocity gradient is dominated by perspec- 
tive rotation, despite the fact that our viewpoint is farther 
from the satellite than the centre of the Galaxy, presumably 
because the satellite is not co-rotating. Velazquez & White 
(1995) made an estimate of 225km s~^ for the tangential 
velocity of the Sagittarius dwarf galaxy by considering only 
the effect of perspective rotation and this rou ghly agrees 



with the observed proper motion of 250km s ^ ( [bata et al 



1997) 



Rotation has been detected in several observational sur- 
veys of Galactic globular clusters. Merritt, Meylan & Mayor 
(1997) analysed the phase-space distribution of 469 stars in 
LoCen and found, after taking perspective rotation into ac- 
count, it to be consistent with solid body rotation out to 
llpc from the centre of the cluster, and decreasing beyond. 
Drukier et al. (1997), in an analysis of the velocities of 230 
stars in the outskirts of M15, found some evidence for ro- 
tation. In both these studies, rotation was found to be less 
important towards the edge of the cluster, suggesting that 
in these cases it is intrinsic and not due to tidal torquing or 
stripping of the bound system. 

4-3.2 Velocity dispersions and tidal heating 

In most panels of Figures |l|and|l| there is good agreement 
between (Tioa calculated with and without the unbound stars. 
However, the observed results can be confusing when looking 
along the orbit (y-axis view panels of Figure ^8|), as there 
is significant unbound material along the line of sight from 
the trailing and leading streamers. This problem may not be 
so severe in reality if the sample can be selected to exclude 
stars beyond a few tidal radii from the cluster (in the figure, 
all stars at the projected separation were included). 

An observer would also see an apparently enhanced ve- 
locity dispersion when looking at the outskirts of the cluster 
from the centre of the Galaxy. In particular, note that the 
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Figure 20. Mass function index a: as a function of time for all 
stars bcund to Model (0a,p3) (solid line), and as measured at 
various joints in the cluster. 
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Figure 21. Initial (solid lines) and final mass and luminosity 
functions for Model (0a:p3) simulated with (dashed lines) and 
without (doti ed lines) diffusion effects. 



dispersion in tiie unbound material roughly corresponds to 
the dispersion of the bound material within the point where 
stripping occurs. This suggests that satellites that are being 
more violently stripped will have larger dispersions in their 
debris trails, though not in excess of the maximum disper- 
sion of the bound material. The Sagittarius dwarf spheroidal 
galaxy is an example of a system where it can be plausibly 
argued that this behaviour is observed: its highly distorted 
surface density contours suggest that it is likely to be sur- 
rounded by a cloud of unbound stars and its velocity disper- 
sion i s roughly constan t along the entire length of its major 
axis (Ibata et al. 1997). 

In their analysis of the velocities in the outskirts of M15, 
Drukier et al. (1997) found that the dispersion of the stars 
decreased to a minimum at 7 arcmin, increasing again be- 
yond this radius. They suggest that this deviation from the 
behaviour expected for an isolated cluster could be due to 
tidal heating, and our simulations confirm this interpreta- 
tion. They comment that the radial position of this mini- 
mum is much smaller than the tidal radius found by G95 
by fitting King Models to star count profiles. This is also 
consistent with our analyses - the star counts and velocity 
analyses become contaminated by unbound stars well within 
the outermost radius of the bound system. Hence the mini- 
mum in the velocity dispersion profile is a good indicator of 
where this contamination becomes important, but does not 
necessarily correspond to the edge of the system. 



4.4 Mass functions 

Jf-.Ji-.l Measuring present day mass functions 
In Figures ^ and 21 we mimic the observational a nalyses 



et al. 1998) and M55 ( ^aggia et al. 1997| ). The solid line in 
Figure ^ shows the evolution of the mass function index 
X for Model (0a,p3). The points indicate what x would be 
measured to be at different points along the orbit and at 
different positions in the cluster if viewed from the centre of 
the Galaxy - within the core radius (filled squares), between 
the core and half-light radii (open squares) and between the 
half-light and tidal radii (stars). Each of these points is cal- 
culated using several thousand stars. In Figure 21 we plot 
the mass and luminosity functions at the beginning (solid 
lines), and the end (dotted lines) of the simulation. The 
dashed lines show the final analysis for the same Model but 
for the simulation that included relaxation effects. 



4.4-^ Placing limits on the initial mass function of a 
cluster 

Suppose we have estimated the fractional mass-loss rate 
from a cluster using its observed population of extra-tidal 
stars. If we also know the PDMF, both globally and lo- 
cally in the outskirts of the cluster, we can place some lim- 
its on its IMF. We test this idea with our simulations by 
using the methods outlined in §4.1 to find (df /dt)\. Figure 
|2^ shows the final mass function (dA^/dM)flnai at the end 
of the five labelled simulations (dashed line) and the mass 
function for the 100 stars at the projected edge of each sys- 
tem [dN / dM)cd^c (dotted line), where each mass function 
has been scaled by the known IMF {{dN/dM)o, represented 
by the horizontal solid line). We then simply estimate the 
IMF (solid squares) to be, 



recently applied to the globular clusters NGC 1261 (Zoccali 



zdATx ^(dN^ (df^ (dN^ 

KdMj o,cst \dMj fl„al ""\dtJl\ dM ) edge 



(20) 
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5.1 Destruction rate of galactic satellites 



o 



0.5 



S 



o 



■ fm ■ ■ ■ ■ 

M- _ m -'-' ■ ! ■ 

■ f ^ fl s ■ 

(l,pi) 



(0a,p3) 



-0.3 



A J! 



-0.4 -0.6 
LOG(M) (MJ 



-0,1 



-1 



Figure 22. Global PDMF (dashed lines) and local PDMF (dot- 
ted lines) for outer 100 stars, for the labelled models. The mass 
functions have been normalized to the IMF of each model (rep- 
resented by solid horizontal lines) and multiplied by arbitrary 
constants. The solid squares show the IMF reconstructed from 
the estimated mass loss rate, global and local mass functions. 
The open squares show the result of the same reconstruction, but 
using the outermost 1000 stars to find the local PDMF. 



where T^h-n is the duration of the simulation. (The open 
squares show the result if the PDMF of the outermost 1000 
stars is used.) 

Several assumptions make this estimate uncertain: (i) 
{df/dt)i is taken from star counts, implicitly assuming that 
the mass function of the extra-tidal material is the same as 
the global mass function; (ii) we assume the mass-function 
of the stripped material to be constant with time; (iii) we 
neglect relaxation effects; and (iv) we assume that the satel- 
lite fills its tidal radius, and hence loses mass at a roughly 
constant rate due to either evaporation or tidal effects. The 
first approximation can be addressed if the luminosity func- 
tion of the cluster is known as a function of radius, as is 
the case in our simulations - Figure demonstrates that 
we can make reasonable predictions for the simulations even 
without this level of refinement; the second two assumptions 
will both lead to an underestimate of the evolution of IMF; 
and the last approximation will place some limits on how 
far back we might be able to 'integrate' the differential mass 
loss. 

Despite these uncertainties, this method provides a new 
approach - more directly based on observations rather than 
using complex dynamical models - to exploring the question 
of whether the IMF in globulars was universal, or environ- 
ment dependent. 



DISCUSSION: APPLICATIONS TO 
OBSERVED SYSTEMS 



5.1.1 Dwarf spheroidal galaxies 

IH95 studied the morphology of eight of the dwarf spheroidal 
satellites of the Milky Way (all the ones currently known 
with the exception of Sagittarius) determined from star 
counts made using the APM facility at Cambridge. Extra- 
tidal stars are clearly seen in the number count profiles for 
six of these satellites, but in no case do they clearly follow 
Sxt ~ r~^. Hence, we use equation (p^, to find an upper 
limit for the mass loss rate. Using the data from Table 3 
of IH95 (kindly made available to us by M. Irwin), we sub- 
tract off the mean background stellar density from each of 
their bins and take rbreak to be the smaller of the point at 
which the slope of the star count profile changes (by visual 
inspection of their figure 2) and the last measured surface 
density above the background. The orbital time period Torb 
is calculated for a circular orbit in a logarithmic potential, 
with circular velocity Vc = 200km s"'^ 



2-kRgc 



1.5 



( Rgc 
\ 50kpc 



Gyrs, 



(21) 



with Rgc taken to be the current distance of each satel- 
lite from the Galactic centre. In fact, if the time period of 
orbits is assumed to be independent of angular momentum 
(shown to be approximately true in Johnston 1998), then 
a satellite at Rgc with velocity v in the range < t; < 
8110/2 will have an orbital time period Torb in the range 
Tcirc(i?Gc)/2 < Torb < 2rcirc(i?Gc), SO equatiou (|l]) is ex- 
pected to be good to within a factor of 2. Several of the 
dwarf spheroidal satellites do have proper motion measure- 
ments - with large errors. These could in principle be used 
to determine a specific time period, but since the mass loss 
estimate itself is expected to be uncertain to within a fac- 
tor of 2, the approximation in equation (ji^) is deemed to 
be sufficiently accurate. For those with proper motions we 
calculate the angle 8 between our line of sight and the plane 
perpendicular to the velocity vector of the satellite in the 
Galactic rest frame. Otherwise, we assume the orbit is cir- 
cular and calculate as the angle between our line of sight 
and the Galactocentric radius vector of the satellite. 

All these quantities are shown in Table ^. The estimates 
{df /dt)2 range from a few percent (Fornax and Leo I) up to 
more than 30 percent in the next Gyr (Ursa Minor and Ca- 
rina). They do not correlate well with previous calculations 
that attempted to determine the robustness of each object 
either from the ratio of the expected to observed tidal radius 
(where the former was calculated given a dynamical estimate 
of the satellite's mass), or of the external tidal field to in- 
ternal field (see IH95). Nevertheless, these mass loss rates 
clearly indicate that the satellite system of the Milky Way 
could easily be diminished by several members in the next 
10 Gyrs, which in turn suggests that there may have been 
several more satellites of the Milky Way in the past. 

If these dSph do indeed have such large mass loss rates, 
why have we not detected tidal streamers? In the cases of 
Ursa Minor and Sculptor this question was addressed by 
Johnston (1998), using a semi-analytic technique. She found 
that if each was losing mass at the rate of 10 percent per 
Gyr the local number count densities along the streamers 
(i.e. not averaged over annuli centred in the cluster) would 
not exceed 1 percent of the background star counts predicted 
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Table 4. Mass loss rates for dwarf spheroidals. 



Name 




'^orb 


e 




(df/dth 




(kpc) 


(Gyr) 


(degrees) 


(Gyr-i) 


Carina 


8.66E+01 


2.60E+00 


2.35EH 


-01 


< 3.31E-01 


Draco 


7.20E+01 


2.16E+00 


4.50EH 


-01 


< 2.18E-01 


Fornax 


1.22E+02 


3.66E+00 


2.94EH 


-01 


< 6.17E-02 


Lcol 


2.02E+02 


6.05E+00 


5.61EH 


-01 


< 5.92E-02 


LeoII 


2.10E+02 


6.29E+00 


3.66EH 


-01 


< 1.44E-01 


Sculptor 


7.22E+01 


2.16E+00 


1.64EH 


-01 


< 2.84E-01 


Sextans 


8.60E+01 


2.58E+00 


4.36E+01 


< 2.61E-01 


Ursa Minor 


6.59E+01 


1.98E+00 


1.30EH 


-01 


< 3.22E-01 



Columns: (1) name; (2) Galactocentric distance; (3) time period of circular orbit at that distance; (4) angle bet ween line of 



s ight and plane perpendicular to satellite 's velocity, calcula ted from proper motion measurements for S culptor ( Schweitzer 
Cudworth, Majewski & Suntzeff (1995)) and Ursa Minor (Schweitzer, Cudworth & Majewski (1998)); (5) mass loss rate 
estimate from equation (|19|). 



by the Bahcall-Soneira (Bahcall & Soneira 198C) model of 
the Milky Way. Hence, even increasing these mass loss rates 
to the tabulated values would not make the streamers strik- 
ing features in the sky. However, the estimates bolster the 
notion that these streamers might be discovered and traced 
over large angular extents using integrated star counts (such 
as the method of Great Circle Cell Counts proposed by John- 
ston, Hernquist & Bolte 1996) or with color and velocity 
information to distinguish them from the background. 



5.1.2 Globular clusters 

Table ^ presents the results of mass loss estimation for the 
12 globular clusters analysed by G95 (using their tables 3- 
14, kindly made available to us by C. Grillmair). Four of 
these clusters have well-defined tidal tails with slopes 7 ~ 1 
(given in column 5), and in these cases both {df/dt)i and 
{df/dt)2 are calculated. In the other cases, only the latter 
estimate is made as an upper limit on the mass loss rate. The 
estimated limits for the mass loss rates range from a few to 
over 100 percent in the next Gyr, again implying that the 
Galaxy's globular cluster system will evolve substantially in 
the next Hubble time. This provides observational support 
for the many purely theoretical studies that ha ve reached 
the same conclusion using semi-analytic models ( Aguilar et 
al. 19881; ICnedin fc Ostriker 1997|; Murali fc Weinberg 1997f 



Vesperini 1997[ [Capriotti fc Hawley 199e| ). 

The last column shows the destruction rates predicted 
by Gnedin & Ostriker (1997). There is no striking corre- 
spondence between our 'observational' results and the purely 
'theoretical' values. However, the upper limit we calculate is 
only in direct contradiction with the theoretical calculation 
in the cases of NGC 288 and NGC 1904. In the former case 
our viewing angle is favourable for making an accurate es- 
timate, but in the latter the value of 9 is sufficiently large 
that we might expect the result to be confused by the de- 
bris geometry. In general, our mass loss rates are higher, but 
Gnedin & Ostriker (1997) themselves point out that their 
destruction rates should be taken as a lower limit since the 
rates could be increased with the inclusion of other effects 
ignored in their study (such as a mass-spectrum). 



5.1.3 Future prospects 

The results in Tables ^ and |^ should be treated with some 
caution as in most cases the observed densities of extra- 
tidal stars do not follow E(r) oc r~^ predicted by the simple 
model which we use to interpret the data. However, the orig- 
inal profiles were made either using star counts directly, or 
with additional photometry to subtract off some of the back- 
ground. Since we expect the tidal debris to typically have 
velocities within ±10kms~^ of the satellite, a study of radial 
velocities around a cluster has the potential of refining these 
measurements considerably. Further in the future, the astro- 
metric satellites SIM (the Space Interferometry Mission) and 
GAIA (the Galactic Astrometric Imaging Satellite) promise 
proper motion measurements with a few jiarcsec accuracy 
(or tangential velocities to better than 1 km s~^ out to tens 
of kpc), and a similar advance in identifying debris. 

We have also, for the sake of simplicity, restricted our 
discussion to annularly averaged surface densities. Clearly, 
more information is contained in two-dimensional surface 



density maps (395). However, these will be more sensitive 
to the orbital phase and the mass loss history of the satel- 
lite and would require more detailed analytic modeling to 
interpret. 



5.2 Determining the IMF from observations of 
the PDMF 

The method proposed in §4.5 uses observations of a cluster's 
PDMF and extra-tidal stars to place some limits on the IMF. 
This is clearly a powerful tool for making progress towards 
understanding whether the IMF is in any sense 'universal'. 

Of course it is non-trivial to find the global PDMF, the 
local mass function in the exterior, and to detect extra-tidal 
stars around a cluster. However, there are currently two ex- 
amples in the literature where this has already been done - 
M15 and M55 (see G95; Piotto, Cool & King 1997; Zaggia 
et al. 1997). In the case of M15, Piotto, Cool & King (1997) 
argue that mass-segregation is important only in the centre 
of the cluster and that their local sample is a fair represen- 
tation of the global PDMF. Unfortunately, in the absence of 
mass segregation our method would find no evolution of the 
IMF since it does not model differential mass-loss due to re- 
laxation effects, but only steady stripping of the most weakly 
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Table 5. Mass loss rates for globular clusters. 



NGC 


Rgc 


^orb 


8 


7 


idf/dt)i 


{df/dt)2 


G&O 




(kpc) 


(Gyr) 


(degrees) 




(Gyr-i) 


(Gyr-i) 


(Gyr-i) 


288 


1.14E+01 


3.43E-01 


4.19E+00 


-5.21E-01 


5.26E-02 


1.26E-02 


l.lOE-01 


362 


9.04E+00 


2.71E-01 


1.40E+01 


-1.58E-01 


6.22E-01 


5.77E-01 


3.54E-02 


1904 


1.81E+01 


5.44E-01 


8.66E+01 






< 8.51E-03 


3.54E-02 


2808 


1.07E+01 


3.20E-01 


2.33E+01 






< 4.37E-02 


1.61E-02 


3201 


8.85E+00 


2.65E-01 


4.95E+01 






< 4.50E-01 


3.45E-02 


4590 


9.94E+00 


2.98E-01 


5.17E-01 






< 1.28E+00 


8.22E-03 


5824 


2.60E+01 


7.81E-01 


7.94E+01 


-1.33E+00 


6.46E-02 


1.18E-01 


3.06E-03 


6864 


1.17E+01 


3.52E-01 


8.60E+01 






< 4.91E-02 


1.89E-02 


6934 


1.17E+01 


3.52E-01 


4.33E+01 






< 4.18E-01 


2.89E-02 


6981 


1.22E+01 


3.65E-01 


5.84E+01 






< 2.96E-01 


1.76E-02 


7078 


l.OlE+01 


3.02E-01 


5.06E+01 






< 7.00E-02 


2.17E-02 


7089 


l.OlE+01 


3.02E-01 


1.91E+01 


-1.62E+00 


1.29E-01 


3.13E-01 


5.76-03 



Columns: (1) name; (2) Galactocentric distance; (3) time period of circular orbit at that distance; (4) angle between line 
of sight and plane perpendicular to satellite's velocity, calculated from proper motions in Dauphole et al. 1996; (5) slope 
of extra-tidal star surface density profile; (6) mass loss rate estimate from equation (p^); (7) mass loss rate estimate from 
equation (|l9[) (8) mass loss rate estimate from Gnedin & Ostriker (1997). 



bound stars. (Since Piotto et al. 1997 make the same state- 
ment about NGC 6397, which has a very different x from 
M15, this might imply that the differences can only be due 
to relaxation effects or an intrinsically different IMF - it is 
unclear whether this is a robust argument, or whether a low 
level of mass segregation farther out in the cluster might be 
sufficient to account for the differences through tidal shock- 
ing.) In the case of M55, Zaggia, Piotto & Capaccioli (1997) 
find some evidence for extra-tidal stars, but do not give the 
density of this material since there are several other effects 
that it could be attributed to. Despite these problems, these 
studies show that it is currently feasible to design future 
observations that can address these issues. 
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6 CONCLUSIONS 

We summarise our main conclusions as follows: 

(i) Observations of extra-tidal stars and enhanced veloc- 
ity dispersions in the outskirts of clusters can be attributed 
to tidally stripped material. The star counts become domi- 
nated by unbound stars at the point where the slope of the 
surface density profile changes or the dispersion reaches a 
minimum. This should not be identified with the tidal ra- 
dius of the cluster since the edge of the bound population 
can still lie significantly beyond this radius. 

(ii) The mass loss rate from a Galactic satellite can be 
directly estimated from the population of extra-tidal stars 
within a few of its tidal radii, the orbital time period and 
our line of sight. 

(iii) Using current observations we calculated the mass- 
loss rate from dwarf spheroidal satellites and globular clus- 
ters and found that both systems will undergo significant 
evolution in the next Hubble time. Future observations 
should be able to place much stronger limits on these de- 
struction rates. 

(iv) Mass-loss estimate together with measurements of 
the local and global PDMFs of a cluster can place limits 
on its IMF. Our review of the literature suggests that it is 
currently observationally feasible to carry out this program. 
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